# create graph for case with interconnection
function create_graph(Rk::Vector{Vector{Int}},ds::DiscreteSpace, Rik::Vector{Vector{Int}},dsi::DiscreteSpace, Ud::Vector{Interval{Float64}}, fxxu::Function)
    dim = ds.dim
    dX = ds.dX

    dimi = dsi.dim
    dXi = dsi.dX

    rtree = create_rtree(Rk,ds)
    G = SimpleDiGraph(length(Rk))
    dd = group_missing_states(Rik)

    @inbounds for (i,bxi) in enumerate(Rk)
        bx = [dX[j][bxi[j]] for j in 1:dim]

        # check if key in dictionary
        if bxi[1] in keys(dd)
            xi =[hull(dXi[1][dd[bxi[1]]])]
            bxy = fxxu(bx,xi,Ud)
            yr = Rect(Tuple(bxy[j].lo for j in 1:dim),Tuple(bxy[j].hi for j in 1:dim))
            bxp = intersects_with(rtree, yr)

            @inbounds for bxpj in bxp
                if i != bxpj.val
                    add_edge!(G,i,bxpj.val)
                end
            end
        end
    end
    return G
end

# create graph for case with interconnection for 4D case
function create_graph_4d(Rk::Vector{Vector{Int}},ds::DiscreteSpace, Rik::Vector{Vector{Int}},dsi::DiscreteSpace, Ud::Vector{Interval{Float64}}, fxxu::Function)
    dim = ds.dim
    dX = ds.dX

    dimi = dsi.dim
    dXi = dsi.dX

    rtree = create_rtree(Rk,ds)
    G = SimpleDiGraph(length(Rk))
    dd = group_missing_states_4d(Rik)

    @inbounds for (i,bxi) in enumerate(Rk)
        bx = [dX[j][bxi[j]] for j in 1:dim]
        # check if key in dictionary
        if (bxi[1],bxi[2]) in keys(dd)
            xii = dd[(bxi[1],bxi[2])]
            xi1 = hull(dXi[1][unique!([i for (i,j) in xii])])
            xi2 =hull(dXi[2][unique!([j for (i,j) in xii])])
            xi = [xi1,xi2]
            bxy = fxxu(bx,xi,Ud)
            yr = Rect(Tuple(bxy[j].lo for j in 1:dim),Tuple(bxy[j].hi for j in 1:dim))
            bxp = intersects_with(rtree, yr)

            @inbounds for bxpj in bxp
                if i != bxpj.val
                    add_edge!(G,i,bxpj.val)
                end
            end
        end
    end
    return G
end

# estimate missing state information and store in dictionary
function group_missing_states(Rk::Vector{Vector{Int}})  
    dd = Dict{Int,Vector{Int}}()
    for (i,j) in Rk
        if j in keys(dd)
            push!(dd[j],i)
        else
            dd[j] = [i]
        end
    end
    return dd
end

# estimate missing state information and store in dictionary for 4D case
function group_missing_states_4d(Rk::Vector{Vector{Int}})  
    dd = Dict{Tuple{Int,Int},Vector{Tuple{Int,Int}}}()
    for (i,j,k,l) in Rk
        if (k,l) in keys(dd)
            push!(dd[(k,l)],(i,j))
        else
            dd[(k,l)] = [(i,j)]
        end
    end
    return dd
end

# find cells to test
function select_test_cells(nlc, lc, G, d, Rkall)

    # lc = setdiff(collect(1:d*d),nlc)
    inn = all_inneighbors(G,lc)
    inn = unique!(vcat(inn...))
    inn = intersect(nlc,inn)
    Rt = Rkall[inn]
    return Rt
end

# reconstruct set
function reconstruct_set(R1k, R2k, Rtt)
    xi = group_missing_states(R2k)

    # reconstruct set
    Rr = unique!([[i,j1,k] for (i,j1) in R1k for (j2,k) in R2k if j1==j2])
    Rt = unique!([[i,j,k] for (j,k) in Rtt for i in xi[j] ])
    return Rr, intersect(Rr,Rt)
end

# validate set
function validate_set(Rr, Rt, ds, Ud, fxu)
    dim = ds.dim
    

    l = 0
    test = true


    while test
        println(l)
        rtree = create_rtree(Rr,ds)
        toremove = []
        for (i,j,k) in Rt
            bx = [ds.dX[1][i], ds.dX[2][j], ds.dX[3][k]]
            bxy = fxu(bx,Ud)
            yr = Rect(Tuple(bxy[j].lo for j in 1:dim),Tuple(bxy[j].hi for j in 1:dim))
            bxp = intersects_with(rtree, yr)

            if isempty(bxp)
                push!(toremove,[i,j,k])
            end
        end

        if isempty(toremove)
            test = false
        end

        Rt = setdiff(Rt,toremove)
        Rr = setdiff(Rr,toremove)
        l = l + 1
    end

    return Rr
end

# validate set 6D
function validate_set_6d(Rr, Rt, ds, Ud, fxu)
    dim = ds.dim
    

    l = 0
    test = true


    while test
        println(l)
        rtree = create_rtree(Rr,ds)
        toremove = []
        for (i,j,k,l,m,n) in Rt
            bx = [ds.dX[1][i], ds.dX[2][j], ds.dX[3][k], ds.dX[4][l], ds.dX[5][m], ds.dX[6][n]]
            bxy = fxu(bx,Ud)
            yr = Rect(Tuple(bxy[j].lo for j in 1:dim),Tuple(bxy[j].hi for j in 1:dim))
            bxp = intersects_with(rtree, yr)

            if isempty(bxp)
                push!(toremove,[i,j,k,l,m,n])
            end
        end

        if isempty(toremove)
            test = false
        end

        Rt = setdiff(Rt,toremove)
        Rr = setdiff(Rr,toremove)
        l = l + 1
    end

    return Rr
end